Small Sample Inference

Simulations to verify assumptions behind various frequentist hypothesis tests and confidence intervals.

Single sample T test

1 January 2024·Updated: 1 January 2025

Based on: Handbook of Parametric and Nonparametric Statistical Procedures, 5th edition, Chapter 2

using Random
Random.seed!(783376894512)
using Distributions
using Intervals
using Plots
n_sim = 1_000_000 # number of simulations to verify theoretical distributions

True model

n = 4
σₜ = 2.0
μₜ = 3.0
Mₜ = Product([Normal(μₜ, σₜ) for _ in 1:n]) # assumed known to be Normal and i.i.d.
generate_data(Mₜ) = rand(Mₜ)

Hypothesis testing

function test_statistic(y, μ₀)
    n = length(y)
    μ̂ = sum(y) / n
    σ̂² = sum((y .- μ̂) .^ 2) / (n - 1)
    return (μ̂ - μ₀) / √(σ̂² / n)
end

Null Hypothesis is true

μ₀ = 3.0
μₜ == μ₀

simulated_statistics = [test_statistic(generate_data(Mₜ), μ₀) for _ in 1:n_sim]
distribution_under_H₀ = t -> pdf(TDist(n - 1), t)

Code

histogram(simulated_statistics; normalize=true, label="Simulated")
tmin = mean(simulated_statistics) -3*std(simulated_statistics)
tmax = mean(simulated_statistics) +3*std(simulated_statistics)
t_plot = tmin:((tmax - tmin) / 1000):tmax
plot!(distribution_under_H₀, t_plot; lw=5, label="Theoretical")
plot!(;
    title="Distribution of single sample T test statistic, \n when Null Hypothesis is correct",
    xlimits=(tmin, tmax),
    xlabel="t",
    ylabel="pdf",
    grid=false,
    legendfontsize=12,
    tickfontsize=12,
    guidefontsize=16,
)

Null Hypothesis is true, but assumptions don’t hold

Observations are non-normal

M_non_normal =  Product([MixtureModel([Normal(μₜ+sqrt(3.99), 0.1),
                                       Normal(μₜ-sqrt(3.99), 0.1),], [0.5, 0.5]) for _ in 1:n])
all(mean(M_non_normal) .== μₜ)
all(sqrt.(var(M_non_normal)) .≈ σₜ)

simulated_statistics = [test_statistic(generate_data(M_non_normal), μ₀) for _ in 1:n_sim]

Code

histogram(simulated_statistics; normalize=true, label="Simulated")
tmin = mean(simulated_statistics) -3*std(simulated_statistics)
tmax = mean(simulated_statistics) +3*std(simulated_statistics)
t_plot = tmin:((tmax - tmin) / 1000):tmax
plot!(distribution_under_H₀, t_plot; lw=5, label="Theoretical")
plot!(;
    title="Distribution of single sample T test statistic, \n under H₀, but non-normal observations",
    xlimits=(tmin, tmax),
    xlabel="t",
    ylabel="pdf",
    grid=false,
    legendfontsize=12,
    tickfontsize=12,
    guidefontsize=16,
)

Observations are non-independent

M_non_independent = MvNormal([μₜ, μₜ, μₜ, μₜ, μₜ], [σₜ^2 0.5*σₜ^2 0 0 0;
                                                0.5*σₜ^2 σₜ^2 0 0 0;
                                                0 0 σₜ^2 0 0;
                                                0 0 0 σₜ^2 0;
                                                0 0 0 0 σₜ^2])
all(mean(M_non_independent) .== μₜ)

simulated_statistics = [test_statistic(generate_data(M_non_independent), μ₀) for _ in 1:n_sim]

Code

histogram(simulated_statistics; normalize=true, label="Simulated")
tmin = mean(simulated_statistics) -3*std(simulated_statistics)
tmax = mean(simulated_statistics) +3*std(simulated_statistics)
t_plot = tmin:((tmax - tmin) / 1000):tmax
plot!(distribution_under_H₀, t_plot; lw=5, label="Theoretical")
plot!(;
    title="Distribution of single sample T test statistic, \n under H₀, but non-independent observations",
    xlimits=(tmin, tmax),
    xlabel="t",
    ylabel="pdf",
    grid=false,
    legendfontsize=12,
    tickfontsize=12,
    guidefontsize=16,
)

Observations are non-identical

M_non_identical = Product([Normal(μₜ, i*σₜ) for i in 1:n])
all(mean(M_non_identical) .== μₜ)

simulated_statistics = [test_statistic(generate_data(M_non_identical), μ₀) for _ in 1:n_sim]

Code

histogram(simulated_statistics; normalize=true, label="Simulated")
tmin = mean(simulated_statistics) -3*std(simulated_statistics)
tmax = mean(simulated_statistics) +3*std(simulated_statistics)
t_plot = tmin:((tmax - tmin) / 1000):tmax
plot!(distribution_under_H₀, t_plot; lw=5, label="Theoretical")
plot!(;
    title="Distribution of single sample T test statistic, \n under H₀, but non-identical observations",
    xlimits=(tmin, tmax),
    xlabel="t",
    ylabel="pdf",
    grid=false,
    legendfontsize=12,
    tickfontsize=12,
    guidefontsize=16,
)

Null Hypothesis is false

μ₀ = 5.0
μₜ != μ₀

simulated_statistics = [test_statistic(generate_data(Mₜ), μ₀) for _ in 1:n_sim]
distribution_under_H₁ = tnc -> pdf(NoncentralT(n - 1, (μₜ - μ₀) / √(σₜ^2 / n)), tnc)

Code

histogram(simulated_statistics; normalize=true, label="Simulated")
tmin = mean(simulated_statistics) -3*std(simulated_statistics)
tmax = mean(simulated_statistics) +3*std(simulated_statistics)
t_plot = tmin:((tmax - tmin) / 1000):tmax
plot!(distribution_under_H₁, t_plot; lw=5, label="Theoretical")
plot!(;
    title="Distribution of single sample T test statistic, \n when Null Hypothesis is incorrect",
    xlimits=(tmin, tmax),
    xlabel="t",
    ylabel="pdf",
    grid=false,
    legendfontsize=12,
    tickfontsize=12,
    guidefontsize=16,
)

Null Hypothesis is false, but assumptions don’t hold

Same counter examples as above.

Single sample Z test

1 January 2024·Updated: 1 January 2025

Based on: Handbook of Parametric and Nonparametric Statistical Procedures, 5th edition, Chapter 1

using Random
Random.seed!(783376894512)
using Distributions
using Intervals
using Plots
n_sim = 1_000_000 # number of simulations to verify theoretical distributions

True model

n = 3
σₜ = 2.0 # assumed known
μₜ = 3.0
Mₜ = Product([Normal(μₜ, σₜ) for _ in 1:n]) # assumed known to be Normal and i.i.d.
generate_data(Mₜ) = rand(Mₜ)

Hypothesis testing

function test_statistic(y, μ₀, σₜ)
    n = length(y)
    μ̂ = sum(y) / n
    return (μ̂ - μ₀) / (σₜ / √n)
end

Null Hypothesis is true

μ₀ = 3.0
μₜ == μ₀

simulated_statistics = [test_statistic(generate_data(Mₜ), μ₀, σₜ) for _ in 1:n_sim]
distribution_under_H₀ = z -> pdf(Normal(), z)

Code

histogram(simulated_statistics; normalize=true, label="Simulated")
zmin, zmax = minimum(simulated_statistics), maximum(simulated_statistics)
z_plot = zmin:((zmax - zmin) / 100):zmax
plot!(distribution_under_H₀, z_plot; lw=5, label="Theoretical")
plot!(;
    title="Distribution of single sample Z test statistic, \n when Null Hypothesis is correct",
    xlimits=(zmin, zmax),
    xlabel="z",
    ylabel="pdf",
    grid=false,
    legendfontsize=12,
    tickfontsize=12,
    guidefontsize=16,
)

Null Hypothesis is true, but assumptions don’t hold

(Known) variance is wrong

σ_wrong = 3.0
M_wrong_variance = Product([Normal(μₜ, σ_wrong) for _ in 1:n])
σ_wrong != σₜ

simulated_statistics = [test_statistic(generate_data(M_wrong_variance), μ₀, σₜ) for _ in 1:n_sim]

Code

histogram(simulated_statistics; normalize=true, label="Simulated")
zmin, zmax = minimum(simulated_statistics), maximum(simulated_statistics)
z_plot = zmin:((zmax - zmin) / 100):zmax
plot!(distribution_under_H₀, z_plot; lw=5, label="Theoretical")
plot!(;
    title="Distribution of single sample Z test statistic, \n under H₀, but variance is wrong",
    xlimits=(zmin, zmax),
    xlabel="z",
    ylabel="pdf",
    grid=false,
    legendfontsize=12,
    tickfontsize=12,
    guidefontsize=16,
)

Observations are non-normal

M_non_normal =  Product([MixtureModel([Normal(μₜ+sqrt(3.99), 0.1),
                                       Normal(μₜ-sqrt(3.99), 0.1),], [0.5, 0.5]) for _ in 1:n])
all(mean(M_non_normal) .== μₜ)
all(sqrt.(var(M_non_normal)) .≈ σₜ)

simulated_statistics = [test_statistic(generate_data(M_non_normal), μ₀, σₜ) for _ in 1:n_sim]

Code

histogram(simulated_statistics; normalize=true, label="Simulated")
zmin, zmax = minimum(simulated_statistics), maximum(simulated_statistics)
z_plot = zmin:((zmax - zmin) / 100):zmax
plot!(distribution_under_H₀, z_plot; lw=5, label="Theoretical")
plot!(;
    title="Distribution of single sample Z test statistic, \n under H₀, but non-normal observations",
    xlimits=(zmin, zmax),
    xlabel="z",
    ylabel="pdf",
    grid=false,
    legendfontsize=12,
    tickfontsize=12,
    guidefontsize=16,
)

Observations are non-independent

M_non_independent = MvNormal([μₜ, μₜ, μₜ], [σₜ^2 0.5*σₜ^2 0; 0.5*σₜ^2 σₜ^2 0; 0 0 σₜ^2])
all(mean(M_non_independent) .== μₜ)
all(sqrt.(var(M_non_independent )) .≈ σₜ)

simulated_statistics = [test_statistic(generate_data(M_non_independent), μ₀, σₜ) for _ in 1:n_sim]

Code

histogram(simulated_statistics; normalize=true, label="Simulated")
zmin, zmax = minimum(simulated_statistics), maximum(simulated_statistics)
z_plot = zmin:((zmax - zmin) / 100):zmax
plot!(distribution_under_H₀, z_plot; lw=5, label="Theoretical")
plot!(;
    title="Distribution of single sample Z test statistic, \n under H₀, but non-independent observations",
    xlimits=(zmin, zmax),
    xlabel="z",
    ylabel="pdf",
    grid=false,
    legendfontsize=12,
    tickfontsize=12,
    guidefontsize=16,
)

Observations are non-identical

There is no way of violating the identical distribution assumption, without also violating another assumption already present in the previous 3 counterexamples.

single sample χ² test for population variance

1 January 2024·Updated: 1 January 2025

Based on: Handbook of Parametric and Nonparametric Statistical Procedures, 5th edition, Chapter 3

using Random
Random.seed!(783376894512)
using Distributions
using Intervals
using Plots
n_sim = 1_000_000 # number of simulations to verify theoretical distributions

True model

n = 4
σₜ = 2.0
μₜ = 3.0
Mₜ = Product([Normal(μₜ, σₜ) for _ in 1:n]) # assumed known to be Normal and i.i.d.
generate_data(Mₜ) = rand(Mₜ)

Hypothesis testing

function test_statistic(y, σ₀)
    n = length(y)
    μ̂ = sum(y) / n
    σ̂² = sum((y .- μ̂) .^ 2) / (n - 1)
    return σ̂² / (σ₀^2 / (n - 1))
end

Null Hypothesis is true

σ₀ = 2.0
σₜ == σ₀

simulated_statistics = [test_statistic(generate_data(Mₜ), σ₀) for _ in 1:n_sim]
distribution_under_H₀ = χ² -> pdf(Chisq(n - 1), χ²)

Code

histogram(simulated_statistics; normalize=true, label="Simulated")
χmin = mean(simulated_statistics) -3*std(simulated_statistics)
χmax = mean(simulated_statistics) +3*std(simulated_statistics)
χ_plot = χmin:((χmax - χmin) / 1000):χmax
plot!(distribution_under_H₀, χ_plot; lw=5, label="Theoretical")
plot!(;
    title="Distribution of Χ² test statistic for a pop. var., \n when Null Hypothesis is correct",
    xlimits=(χmin, χmax),
    xlabel="χ",
    ylabel="pdf",
    grid=false,
    legendfontsize=12,
    tickfontsize=12,
    guidefontsize=16,
)

Null Hypothesis is true, but assumptions don’t hold

Observations are non-normal

M_non_normal =  Product([MixtureModel([Normal(μₜ+sqrt(3.99), 0.1),
                                       Normal(μₜ-sqrt(3.99), 0.1),], [0.5, 0.5]) for _ in 1:n])
all(mean(M_non_normal) .== μₜ)
all(sqrt.(var(M_non_normal)) .≈ σₜ)

simulated_statistics = [test_statistic(generate_data(M_non_normal), σ₀) for _ in 1:n_sim]

Code

histogram(simulated_statistics; normalize=true, label="Simulated")
χmin = mean(simulated_statistics) -3*std(simulated_statistics)
χmax = mean(simulated_statistics) +3*std(simulated_statistics)
χ_plot = χmin:((χmax - χmin) / 1000):χmax
plot!(distribution_under_H₀, χ_plot; lw=5, label="Theoretical")
plot!(;
    title="Distribution of Χ² test statistic for a pop. var., \n under H₀, but non-normal observations",
    xlimits=(χmin, χmax),
    xlabel="χ",
    ylabel="pdf",
    grid=false,
    legendfontsize=12,
    tickfontsize=12,
    guidefontsize=16,
)

Observations are non-independent

M_non_independent = MvNormal([μₜ, μₜ, μₜ, μₜ, μₜ], [σₜ^2 0.9*σₜ^2 0 0 0;
                                                0.9*σₜ^2 σₜ^2 0 0 0;
                                                0 0 σₜ^2 0 0;
                                                0 0 0 σₜ^2 0;
                                                0 0 0 0 σₜ^2])

all(sqrt.(var(M_non_independent)) .≈ σₜ)

simulated_statistics = [test_statistic(generate_data(M_non_independent), σ₀) for _ in 1:n_sim]

Code

histogram(simulated_statistics; normalize=true, label="Simulated")
χmin = mean(simulated_statistics) -3*std(simulated_statistics)
χmax = mean(simulated_statistics) +3*std(simulated_statistics)
χ_plot = χmin:((χmax - χmin) / 1000):χmax
plot!(distribution_under_H₀, χ_plot; lw=5, label="Theoretical")
plot!(;
    title="Distribution of Χ² test statistic for a pop. var., \n under H₀, but non-independent observations",
    xlimits=(χmin, χmax),
    xlabel="χ",
    ylabel="pdf",
    grid=false,
    legendfontsize=12,
    tickfontsize=12,
    guidefontsize=16,
)

Observations are non-identical

The only way to create non-identical observations, while keeping the other assumptions, is to vary the mean.

Two independent samples T test

1 January 2024·Updated: 1 January 2025

Based on: Handbook of Parametric and Nonparametric Statistical Procedures, 5th edition, Chapter 11

using Random
Random.seed!(783376894512)
using Distributions
using Intervals
using Plots
using LinearAlgebra
n_sim = 1_000_000 # number of simulations to verify theoretical distributions

True model

n1 = 4
n2 = 5
σₜ = 2.0
μₜ1 = 3.0
μₜ2 = 4.0
Mₜ1 = Product([Normal(μₜ1, σₜ) for _ in 1:n1]) # assumed known to be Normal and i.i.d.
Mₜ2 = Product([Normal(μₜ2, σₜ) for _ in 1:n2]) # assumed known to be Normal and i.i.d.
generate_data(M1, M2) = (rand(M1), rand(M2))

Hypothesis testing

function test_statistic(y1, y2, Δμ₀)
    n1 = length(y1)
    n2 = length(y2)
    μ̂1 = sum(y1) / n1
    μ̂2 = sum(y2) / n2
    σ̂²1 = sum((y1 .- μ̂1) .^ 2) / (n1 - 1)
    σ̂²2 = sum((y2 .- μ̂2) .^ 2) / (n2 - 1)
    σ̂²p = ((n1 - 1) * σ̂²1 + (n2 - 1) * σ̂²2) / (n1 + n2 - 2)
    return (μ̂1 - μ̂2 - Δμ₀) / (sqrt(σ̂²p) * sqrt(1 / n1 + 1 / n2))
end

Null Hypothesis is true

Δμ₀ = -1.0
Δμ₀ == μₜ1 - μₜ2

all(mean(Mₜ1) .== μₜ1)
all(mean(Mₜ2) .== μₜ2)
all(sqrt.(var((Mₜ1))) .== σₜ)
all(sqrt.(var((Mₜ2))) .== σₜ)

simulated_statistics = [test_statistic(generate_data(Mₜ1, Mₜ2)..., Δμ₀) for _ in 1:n_sim]
distribution_under_H₀ = t -> pdf(TDist(n1 + n2 - 2), t)

Code

histogram(simulated_statistics; normalize=true, label="Simulated")
tmin = mean(simulated_statistics) - 3 * std(simulated_statistics)
tmax = mean(simulated_statistics) + 3 * std(simulated_statistics)
t_plot = tmin:((tmax - tmin) / 1000):tmax
plot!(distribution_under_H₀, t_plot; lw=5, label="Theoretical")
plot!(;
    title="Distribution of two independent samples T test statistic, \n when Null Hypothesis is correct",
    xlimits=(tmin, tmax),
    xlabel="t",
    ylabel="pdf",
    grid=false,
    legendfontsize=12,
    tickfontsize=12,
    guidefontsize=16,
)

Null Hypothesis is true, but assumptions don’t hold

Observations are non-normal

M1_non_normal = Product([
    MixtureModel(
        [Normal(μₜ1 + sqrt(3.99), 0.1), Normal(μₜ1 - sqrt(3.99), 0.1)], [0.5, 0.5]
    ) for _ in 1:n1
])
M2_non_normal = Product([
    MixtureModel(
        [Normal(μₜ2 + sqrt(3.99), 0.1), Normal(μₜ2 - sqrt(3.99), 0.1)], [0.5, 0.5]
    ) for _ in 1:n2
])
all(mean(M1_non_normal) .== μₜ1)
all(mean(M2_non_normal) .== μₜ2)
all(sqrt.(var(M1_non_normal)) .== σₜ)
all(sqrt.(var(M2_non_normal)) .== σₜ)
simulated_statistics = [
    test_statistic(generate_data(M1_non_normal, M2_non_normal)..., Δμ₀) for _ in 1:n_sim
]

Code

histogram(simulated_statistics; normalize=true, label="Simulated")
tmin = mean(simulated_statistics) - 3 * std(simulated_statistics)
tmax = mean(simulated_statistics) + 3 * std(simulated_statistics)
t_plot = tmin:((tmax - tmin) / 1000):tmax
plot!(distribution_under_H₀, t_plot; lw=5, label="Theoretical")
plot!(;
    title="Distribution of two independent samples T test statistic, \n under H₀, but non-normal observations",
    xlimits=(tmin, tmax),
    xlabel="t",
    ylabel="pdf",
    grid=false,
    legendfontsize=12,
    tickfontsize=12,
    guidefontsize=16,
)

Observations are non-independent

M1_non_independent = MvNormal(
    [μₜ1, μₜ1, μₜ1, μₜ1, μₜ1],
    [
        σₜ^2 0.5*σₜ^2 0 0 0
        0.5*σₜ^2 σₜ^2 0 0 0
        0 0 σₜ^2 0 0
        0 0 0 σₜ^2 0
        0 0 0 0 σₜ^2
    ],
)
Σ2 = collect(Diagonal(fill(σₜ^2, n2)))
Σ2[2, 1] = Σ2[1, 2] = 0.5 * σₜ^2
M2_non_independent = MvNormal(fill(μₜ2, n2), Σ2)

all(mean(M1_non_independent) .== μₜ1)
all(mean(M2_non_independent) .== μₜ2)
all(sqrt.(var(M1_non_independent)) .== σₜ)
all(sqrt.(var(M2_non_independent)) .== σₜ)
simulated_statistics = [
    test_statistic(generate_data(M2_non_independent, M2_non_independent)..., Δμ₀) for
    _ in 1:n_sim
]

Code

histogram(simulated_statistics; normalize=true, label="Simulated")
tmin = mean(simulated_statistics) - 3 * std(simulated_statistics)
tmax = mean(simulated_statistics) + 3 * std(simulated_statistics)
t_plot = tmin:((tmax - tmin) / 1000):tmax
plot!(distribution_under_H₀, t_plot; lw=5, label="Theoretical")
plot!(;
    title="Distribution of two independent samples T test statistic, \n under H₀, but non-independent observations",
    xlimits=(tmin, tmax),
    xlabel="t",
    ylabel="pdf",
    grid=false,
    legendfontsize=12,
    tickfontsize=12,
    guidefontsize=16,
)

Observations are non-identical

M1_non_identical = Product([Normal(μₜ1, i * σₜ) for i in 1:n1])
M2_non_identical = Product([Normal(μₜ2, i * σₜ) for i in 1:n1]) # note the n1 here to differentiate from the next counterexample
all(mean(M1_non_identical) .== μₜ1)
all(mean(M2_non_identical) .== μₜ2)
any(sqrt.(var(M1_non_identical)) .== σₜ)
any(sqrt.(var(M2_non_identical)) .== σₜ)

simulated_statistics = [
    test_statistic(generate_data(M1_non_identical, M2_non_identical)..., Δμ₀) for
    _ in 1:n_sim
]

Code

histogram(simulated_statistics; normalize=true, label="Simulated")
tmin = mean(simulated_statistics) - 3 * std(simulated_statistics)
tmax = mean(simulated_statistics) + 3 * std(simulated_statistics)
t_plot = tmin:((tmax - tmin) / 1000):tmax
plot!(distribution_under_H₀, t_plot; lw=5, label="Theoretical")
plot!(;
    title="Distribution of two independent samples T test statistic, \n under H₀, but non-identical observations",
    xlimits=(tmin, tmax),
    xlabel="t",
    ylabel="pdf",
    grid=false,
    legendfontsize=12,
    tickfontsize=12,
    guidefontsize=16,
)

Two samples have a different variance

M2_different_variance = Product([Normal(μₜ1, 2 * σₜ) for _ in 1:n2])
all(mean(M2_different_variance) .== μₜ1)
all(sqrt.(var(M2_different_variance)) .!= σₜ)

simulated_statistics = [
    test_statistic(generate_data(Mₜ1, M2_different_variance)..., Δμ₀) for _ in 1:n_sim
]

Code

histogram(simulated_statistics; normalize=true, label="Simulated")
tmin = mean(simulated_statistics) - 3 * std(simulated_statistics)
tmax = mean(simulated_statistics) + 3 * std(simulated_statistics)
t_plot = tmin:((tmax - tmin) / 1000):tmax
plot!(distribution_under_H₀, t_plot; lw=5, label="Theoretical")
plot!(;
    title="Distribution of two independent samples T test statistic, \n under H₀, but different sample variances",
    xlimits=(tmin, tmax),
    xlabel="t",
    ylabel="pdf",
    grid=false,
    legendfontsize=12,
    tickfontsize=12,
    guidefontsize=16,
)

Null Hypothesis is false

Δμ₀ = 2.0
Δμ₀ != μₜ1 - μₜ2

simulated_statistics = [test_statistic(generate_data(Mₜ1, Mₜ2)..., Δμ₀) for _ in 1:n_sim]
distribution_under_H₁ =
    tnc -> pdf(NoncentralT(n1 + n2 - 2, (μₜ1 - μₜ2 - Δμ₀) / √(σₜ^2 / n1 + σₜ^2 / n2)), tnc)

Code

histogram(simulated_statistics; normalize=true, label="Simulated")
tmin = mean(simulated_statistics) - 3 * std(simulated_statistics)
tmax = mean(simulated_statistics) + 3 * std(simulated_statistics)
t_plot = tmin:((tmax - tmin) / 1000):tmax
plot!(distribution_under_H₁, t_plot; lw=5, label="Theoretical")
plot!(;
    title="Distribution of two independent samples T test statistic, \n when Null Hypothesis is incorrect",
    xlimits=(tmin, tmax),
    xlabel="t",
    ylabel="pdf",
    grid=false,
    legendfontsize=12,
    tickfontsize=12,
    guidefontsize=16,
)

Null Hypothesis is false, but assumptions don’t hold

Same counter examples as above.